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^ Tlie problem of stationary heat transport in the Fermi-Pasta-Ulam chain is numerically 

studied showing that the conductivity diverges in the thermodynamic limit. Simulations were 
performed with time-reversible thermostats, both for small and large temperature gradients. 
In the latter case, fluctuations of the heat current are shown to be in agreement with the recent 
conjectures of Gallavotti and Cohen [Phys. Rev. Lett. 74, 2694 (1995)]. 
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^ ' I. INTRODUCTION 

in 

^ , 

, The problem of heat conduction in one-dimensional insulating solids is very old. The celebrated Peierls' 
' theory was successful in describing the low-temperature dependence of the thermal conduction cofRcient, 
C clarifying the basic role of Umklapp scattering induced by nonlincarity. However, for sufficiently high ener- 

gies/temperatures, where the usual picture of weakly interacting phonons is no longer appropriate, one has to 
face a fully nonlinear problem and a complete analytical solution seems hardly feasible. 

Accordingly, in such a regime, the problem has been attacked several times by means of nonequilibrium 
molecular-dynamics (see Ref. [|| and the bibliography therein). This amounts to solving numerically the mi- 
croscopic (classical) equations of motion, with the principal goals of recovering the macroscopic heat-diffusion 
fH law and of measuring the corresponding transport coefficient. However, even this program turned out to be far 
Q from trivial, especially in one spatial dimension, where the interpretation of the numerical analysis appears to 
O be rather controversial. 

^ I The only clear results refer to some ID toy models, characterized by a non-smooth interaction, like the so 
called ding-a-ling model ||^ and its simplified version studied in Ref. Q . The remarkable success of such artificial 
systems in producing a "normal" (i.e. length-independent) heat conductivity is claimed to follow from the strong 
inhibition of coherent soliton propagation (see, however, Ref. |^ for a counterexample). Accordingly, the main 
mechanism at work is the diffusive transport of the energy induced by the underlying chaotic dynamics, which 
thus ensures the validity of the Fourier law. 

As soon as one turns to the more realistic case of smooth interaction potentials, the scenario is definitely more 
ambiguouos. For instance, it is claimed that the conductivity of a diatomic Toda lattice is finite or diverging 
depending on the mass ratio . For the Fermi-Pasta-Ulam (FPU) chain, the results reported in Ref. and 
in this paper do not show any evidence of a normal conductivity. In every case, there is no clear explanation of 
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the observed phenomenology, even in the hmit of smaU apphed gradients, so that it is worth reconsidering the 
whole problem. 

This is not, however, our only motivation. Recently, the approach to non-equilibrium statistical mechanics 
through the introduction of time-reversible thermostats proved to be rather effective in several test problems 
0. In fact, if a system, besides being time-reversible, is "sufficiently chaotic", the tools developed for strictly 
hyperbolic systems allow for some statistical predictions even far from equilibrium Q . Relevant numerical tests 
have been performed on shear-flow |^ and electrical conduction models In the second part of the present 
work, we successfully test these ideas also for the heat conduction problem. 

More specifically, in section II we introduce the model equations, while section III is devoted to a brief 
description of the microscopic quantities of interest. Macroscopic properties of heat conduction are discussed 
in section IV, and the predictions of the "chaotic hypothesis" are tested in section V. The open problems and 
future perspectives are shortly summarized in section V. Finally, the expression for the heat flux measured in 
the simulations is derived in the Appendix. 



II. MODELS WITH TIME-REVERSIBLE THERMOSTATS 



We consider a chain of N anharmonic oscillators, denoting with qi the displacement of the Z-th particle from its 
equilibrium position. The lattice has fixed boundary conditions, go = Qn+i = . In the bulk {I — 2, . . . ,N —1), 
the equations of motion are purely hamiltonian, namely 

qi ^ fi- fi+i , (1) 

where /; — —V'{qi — qi-i) is the force resulting from the interaction potential V (the prime denotes derivative 
with respect to the argument) acting between neighbouring particles. In the following we will consider the FPU 
potential 

Vi^) = 1 +P\ (2) 

with (3 — We anticipate that the results described here in the following do not appear to be peculiar of this 
choice of the potential. A similar scenario has been found for different choices of V . 

The first and the A^-th oscillators interact with two heat reservoirs operating at different temperatures, 
and Tr, respectively (without loss of generality, we assume that Tl > Tr), in order to induce a heat flux through 
the chain. Their equations of motion are 

qi = -Cl qi + fi- h 

qN— ~Cr qN + In ~ In+i ■ (3) 

The action of the thermostats is microscopically modelled by the "thermal" variables (l, Cr, which evolve 
according to the Nose-Hoover dynamical equations 

<.^^(|-l) , (4) 

where the time Q is the thermostat response time. The above prescriptions imply that the kinetic energy of 
the boundary particles fluctuates around the imposed average value, thus simulating a "canonical" dynamics. 
In the limit case — > 0, one has the so-called isokinetic (or Gaussian) thermostat: the kinetic energy is exactly 
conserved and the action of the thermal bath is properly described without the need to introduce a further 
dynamical variable, since C becomes an explicit function of the qs 0. We shall see that there is a price to pay 
for the simplification of the dynamical equations: a larger thermal resistance at the boundaries. 

For all values of O, the equations of motion are left invariant under time-reversal composed with the involution 
X defined as qi —qi, C —( (notice that the ^s are momentum-like quantities . For a discussion of 
time reversibility in such thermostatted models and its implications to nonequilibrium statistical mechanics see 
Ref. [O. 
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III. MICROSCOPIC DEFINITION OF THE THERMODYNAMIC QUANTITIES 



Let us define the local energy density of the chain as h{x,t) — '^ihiS{x — xi), where xi = la + qi is the 
position of the Z-th particle in a lattice of spacing a, while 

»; 2 1 

hr.= ^ +^[V{qi+i-qi) + Viqi-qi-l)] (5) 

is the energy per particle (as usual, all the masses are set equal to 1 so that pi = qi). If local equilibrium holds, 
the definition of kinetic temperature stems from the local version of the "virial theorem" 

dhi\ / dhA 



dpi I 



where (■) denotes the time average over a sufficiently long trajectory. Assuming, as it is indeed the case of our 
model, that Eq. (^) holds, we define T; = (p^). 

The local heat flux j(x, t) — jiS{x — xi) is implicitely defined by the continuity equation 

A(x,t) + div j{x,t) = . (7) 

The very definition of ji is non-trivial in strongly anharmonic systems i.e. at high energies/temperatures. The 
correct way to proceed (see the Appendix and Ref . ) is to perform the spatial Fourier transform of Eq. (0) 
and to expand the result in powers of the wavenumber k. Neglecting all higher-order terms, as it is usually 
done in hydrodynamics, one eventually obtains that the heat flux at the l-th position is given by 

Ji{t) ^ lapi{fi+, + fi) , (8) 

so that pifi+i has the simple interpretation of the flow of potential energy from the l-th to the neighbouring 
particle. In the simulations, we have measured the time average {ji) in the bulk. Notice that the stationarity 
condition implies {pifi+i) = (pifi), as it can be easily derived from the condition (qi) = 0. 

For later convenience wc define the total flux through a chain or a subchain of N particles, as the integral of 
j{x,t), namely 



At) = ^E^'W ■ (9) 

At equilibrium, the instantaneous fluxes fluctuate "randomly" around zero. 



IV. MACROSCOPIC PROPERTIES IN THE NONEQUILIBRIUM STATE 

Here, we describe the outcome of our molecular-dynamics simulations, performed with fixed values of the 
temperatures Tr and T^. With this setting, larger numbers of particle correspond to smaller temperature 
gradients (i.e. small external fields), so that, for N — > oo, equilibrium dynamics is locally approached along 
the chain. In other words, we are considering cases where the usual linear response or Green-Kubo theory is 
eventually applicable. Nevertheless, the problem remains highly nontrivial and deserves, as we will see, much 
attention. 

The numerical simulations have been performed with several values of N up to 4096, monitoring both the 
kinetic temperatures and the heat fluxes. In every case, after a suitably long transient, the system reaches a 
statistically stationary state, where each oscillator is in local equilibrium at a certain kinetic temperature and 
Eq. (^ is well verified. 



A. The temperature field 



Whenever the imposed temperatures are different from one another, the microscopic dynamics corresponds 
to a nonequilibrium macrostate characterized by a nonuniform temperature field along the chain (see Fig. 1). 
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FIG. 1. Temperature profiles with Tl = 152, Tr = 24, 6 = 1, for different lattice lengths. Note that 
a nonlinear profile sets in for the larger A*' values. For A'^ = 66 we report the temperature fields for two 
other values of the response time, namely = 0.1 (dashed line) and O = 10 (dot-dashed line). 

Long averaging reveals that the temperature gradient is rather smooth except at the boundaries, where thermal 
resistance effects may generate large temperature jumps. Such effects turn out to depend on the response time 
of the thermostats and are particularly relevant for small Q's (see third panel in Fig. 1). One can understand 
this phenomenon by realizing that a fast reaction of the thermostats makes the dynamics of the end particles 
qualitatively different from that of the neighbouring ones. As one is interested in measuring the bulk contribution 
to the conductivity, it is important to minimize the boundary effects, i.e. the temperature gaps observed at both 
extrema]^ This implies that 8 should be chosen as large as possible; however, the larger is 8, the longer must 
last the simulations in order to have reliable statistics. We find that the choice 8 = 1 represents a reasonable 
compromise. 

As already observed in Rcf. |0| , the behaviour of the temperature profiles for different numbers of particles is 
well reproduced by the scaling Ansatz, 



Tl = T{l/N) 



(10) 



i.e., the shape of the profile is independent of N, if the chain length is rescaled to 1. This implies that 
the temperature field scales everywhere in the same manner and one can choose equivalently any interval for 
measuring the temperature gradient (provided that the interval is measured in fractions of the total length). 

The temperature profile in the bulk is not exactly linear (except for rather small Ns), as one could expect 
from the stationary solution of the heat equation (namely the Fourier law). The nonlinearity of the profiles 
could at first be interpreted as an indication of a temperature-dependent conductivity, but this is actually not 
the case. In fact, simulations performed with as small temperature-differences as Tl — Tr = 4 and for N = 128, 
still reveal clear deviations from linearity. 



■'"Similar problems arise also with random thermostats; ingeniuous tricks have been worked out to circumvent them Q. 
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FIG. 2. Behaviour of the thermal conductivity, Eq. (|ll|), as a function of the lattice length A*' for 
Tl = 152, Tr = 24. The flux j is computed averaging over one long trajectory (~ lO"), started from 
random initial conditions and after discarding a transient (~ 10^). The inset shows the results of Ref. 



B. Thermal conductivity 

Another relevant macroscopic feature is the onset of a constant, on average, heat flux, namely (j;) = j, for 
every ^ = 2, . . . , — 1 (the flux through the oscillators in contact with the heat reservoirs must account also for 
the energy exchange with each thermostat - see the following subsection). 

The thermal conductivity is defined (to the lowest order in the applied gradient) as the ratio 

Our main result is that j , which vanishes as — s- 00, approaches more slowly than the temperature gradient, 
thus implying a diverging conductivity in the thermodynamic limit. In fact, the simulations do reveal that the 
heat flux scales as j oc N~"-, with a definitely smaller than 1. Accordingly (see Fig. 2), k, diverges as N^~°', 
since the scaling behaviour of the profile implies that dT/dx is proportional to (T^ — Tji)/N . A careful scrutiny 
of the data in Fig. 2 reveals a sort of crossover from a « 0.55, for N < 250 to a « 0.62 for larger numbers of 
particles. 

It is worth mentioning that a similar behaviour of the conductivity is found also by using other thcrmostatting 
schemes, such as stochastic interactions with a gas of Maxwellian particles. For comparison, in the inset of Fig. 2, 
we report the data taken from [Q. Although they refer to different temperatures (Tl = 1500, Tr = 150, if 
expressed with reference to the same value of /3 here employed), the scaling behaviour is approximately the 
same as ours. 

Therefore, we are forced to conclude that, at least for the considered system sizes, Fourier law is not satisfied 
and chaoticity is not sufficient to ensure its validity. 



C. Entropy production 

The energy flux at the chain ends is simply given by 

iL,B^ = -{C,L,B,qi,N) ^ ~{Cl,r)Tl.r , (12) 

where the second equality is obtained from the condition {dC,\ ji/dt) — 0. In the stationary regime, the balance 
between the ingoing and outgoing fluxes implies that Jl — —jn which, in turn, implies that (Cl) and {C_r) 
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must have opposite signs. This rather obvious relationship, stemming from energy conservation, has a somehow 
surprising meaning, if interpreted from the point of view of dynamical equations. In fact, a negative (^l) 
(the flux must be obviously positive in the left-end, which is in contact with the hotter reservoir) means an 
expansion of volumes rather than a dissipation! The apparent anomaly is immediately clarified by noticing that 
the system in the whole is globally dissipative as 7 = {(r + Cl) turns out to be positive in all simulations that 
we have performed. This is also consistent with a theorem recently proven by Ruelle [ p5[ for time-reversible 
systems and expresses the intuitive notion that if the energy is conserved on the average, it is not possible that 
volumes steadily diverge. What is not a priori trivial is that a finite dissipation spontaneously arises as soon as 
a non-equilibrium stationary state sets in ||^,^ . 

Eq. ( [l2|) is also susceptible of an interpretation in terms of entropy production. By subtracting from one 
another the two expressions of the fluxes and noticing that j = Jl = — J_r > 0, one obtains 

(Cl) + {(r) . (13) 

which can be interpreted as a balance relation for the global entropy production. In fact, according to the 
principles of irreversible thermodynamics, the local rate of entropy production a in the bulk is given by 

Upon integrating Eq. (p^, the r.h.s of Eq. (^3|) is obtained, which can thus be interpreted as the global 
production rate of entropy in the bulk. On the other hand, according to general arguments on reversible 
thermostats g, the l.h.s. of Eq. (|l|) can be identified with the entropy production from the heat baths. 
Eq. ( p^ has been accurately verified in a wide range of temperatures. 



V. CHAOTIC HYPOTESIS AND LARGE FLUCTUATIONS 



We now turn our attention to far-from-equilibrium regimes, where the applied gradient is very large, i.e. to 
the case of short chains with large temperature differences at the boundaries. Our aim is to check the validity 
of the fluctuation theorem, recently proposed by Gallavotti and Cohen Q , and carefully tested in the periodic 
Lorentz gas with an electric applied field |l^. The test is rather crucial for at least two reasons: (i) at variance 
with most of the cases considered in literature, here only the boundary particles are thermostatted ||l6|; (ii) it is 
not a priori clear that the chaotic hypothesis, one of the key assumptions for the validity of the theorem, applies 
to the present case. In fact, it is far from obvious that a dynamical system such as an FPU chain, displaying a 
slow convergence to equipartition at low temperatures, is "Anosov-like" in the thermodynamic limit. 

The fluctuation theorem essentially connects the probability of positive to that of negative values of the 
entropy production. Although we recommend the reader to consult |^ for a detailed exposition of the theorem, 
here we summarize the main lines of the proof to give a flavour of the various steps. The starting observation is 
that any two regions in phase space|^, mutually related by the involution X, are characterized by opposite entropy 
values a and —a (the entropy changes sign upon applying the involution I). The second key point is that the 
existence of a Sinai-Ruelle-Bowen measure implies that the probability of one such region R is proportional 
to the product of the expanding multipliers (over a suitable time lag). As a consequence, invariance under 
time-reversal ensures that the probability to observe —a in R is also proportional to the product of the inverse 
of the contracting multipliers in I{R) . This fact implies that the ratio of the probability of observing a to that 
of —a can be reduced to the product of expanding and contracting multipliers, all taken at the point where a is 
actually observed. This product is, in turn, nothing but the volume contraction or, equivalently, the exponential 
of the entropy production. 

From Eq. (ES), one can see that, apart from an irrelevant multiplicative constant (the external- field amplitude), 
the heat flux j is equivalent to the entropy production. In order to make the analogy with the previous study 
even more stringent, we have, however, preferred to consider the global heat flux as defined in Eq. (^ as the 
latter quantity corresponds to a spatial average over all oscillators (let us indeed recall that in many previous 



Technically speaking, we should refer to a sufficiently small element of a Markov partition. 
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FIG. 3. The distribution Pt{z) for different values r. The chain lenght is A*" = 14 and the boundary 
temperatures are Tl ~ 120, Tr — 56 and tq — 50. The value of the average flux is j — 30.67 ± 0.01. 



studies the entropy production is an extensive quantity, proportional to the number of degrees of freedom). 
This change of definition does not affect the symmetry properties of J, so that the fluctuation theorem applies 
(or not) independently of this modification. 

Therefore, our starting point is the finite-time average of the global heat flux (proportional to the entropy 
production) , 

{J)r = - r\j{t')dt' , (15) 



where the average is performed over a sufficiently long trajectory to ensure a good accuracy of the underlying 
Markov approximation (see again Rcf. Q). We then compute the probability distribution Pr of the variable 

- - ^'^^ (16) 



where (J)oo, denoting the "infinite time" average of the global flux, coincides with the stationary average flux 
j introduced in the previous sections. The Gallavotti-Cohcn conjecture then reads as 



Pr{-z) \Tr Tl 

We have performed simulations for N — 14, Tl — 120, Tr — 56, = 1 and several values of r ranging from 5 to 
80. At variance with the results of Ref. ||l^, the bell-shaped distribution P{z) (Fig. 3) is clearly not Gaussian: 
both of its tails approach zero exponentially for large values of |z|, but with different rates. Upon increasing r, 
the central part of Pt{z) approaches more and more a Gaussian shape as one would expect from the increasing 
statistical independence of the data. On the other hand, reliable numerical measures become more and more 
difficult: for instance, already for t = 80, negative values of z are much too rare to permit the evaluation of 
the distribution over a reasonable integration time. Nevertheless, the linear behaviour predicted by Eq. (|l7| ) is 
indeed observed, as shown in Fig. 4, where the numerical points are superimposed to the theoretical prediction. 
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FIG. 4. Test of the fluctuation theorem, for the same parameters as in Fig. 3 and different values of the 
time r. The sohd hnes are the predictions of Eq. (p^). 

VI. CONCLUSIONS AND PERSPECTIVES 

In this paper we have shown that the fluctuation theorem recently proposed by Gallavotti and Cohen |^ is 
successfuUy verified for a chain of FPU oscillators in contact with time-reversible thermostats at the boundaries. 
It is worth stressing that this is a very different context with respect to the ones previously considered in the 
literature In this sense, our results definitely enforce the generality of this approach in describing the 

thermodynamics of nonequilibrium stationary states in terms of the Sinai-Bowen-Ruelle probability distribution. 

Furthermore, we have found that heat conductivity appears to diverge in the thermodynamic limit. This is 
at variance with results obtained in with reference to the so-called ding-a-ling model. The origin of the 
divergence of k in FPU seems to be traceable to the existence of quasi-conserved long-wavelength modes jT^ 
supporting almost ballistic transport along the chain. Such a feature is indeed generic for Hamiltonian models 
with smooth nearest-neighbour coupling like (|l|), but it is clearly absent in models such as ding-a-ling which 
cannot sustain long- wavelength modes. 

We plan to proceed towards a more detailed comprehension of heat conduction by exploring the following 
routes: 

• Dependence of transport properties on the physical settings; for example, by changing the number of 
particles in thermal contact with the reservoirs. 

• Measurement of the thermal conductivity in the framework of Green-Kubo linear response theory with 
equilibrium simulations. 

• Application of the same approaches to different models of nonlinear chains of oscillators, including the 
extension to two and three dimensions. 

In particular, we want to stress the importance of the last point since it is not unlikely that macroscpic 
validity of the Fourier law is generally ensured only in some higher dimension. After all, statistical mechanics 
is full of phenomena that do depend on the dimensionality. 
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APPENDIX A: 

In this Appendix we detemine the expression of the heat flux for a Id chain, where the neighbouring particles 
interact via a generic nonlinear potential, such that the resulting equations of motion arc of the form (|l|). For 
simplicity we assume that the chain is infinite. 

Let us write the energy density and the flux as Fourier integrals 

h{k,t) = j h{x,t)e-'''''dk j{k,t) ^ J j{x,t)e-'''''dk , (Al) 
so that Eq. (0) becomes 

^+*fcj = . (A2) 
We now rewrite the last equation, by splitting the heat flux into two different contributions 

J = , (A3) 

such that, by the explicitely computing the time-derivative of H, it is recognized that 

-^pz/iie-*'^^' , (A4) 

while 

As we want to determine j^^^ to the lowest order in ik, we expand the exponentials as 
end then 

Summing the two terms and recalling the definition of j, it can be realized that (after some further index 
manipulation) 

ji = Pihi + ^{qi - qi-i)ipi +pi-i)fi + ^pi{fi + fi+i) . (A8) 

With similar calculations, it can be shown that the first two terms in the r.h.s. correspond to order-fc'^ terms in 
the Fourier expansion so that they can be neglected and Eq. (H) is finally obtained. 
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